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Abstract-  A  model-driven,  multiscale  medical  image 
segmentation  system  is  presented.  A  tree  representation  is 
calculated  for  the  image,  using  a  modification  of  the  immersion 
algorithm  used  for  watersheds  calculation.  Segmentation  is 
carried  out  by  a  matching  process  between  the  obtained  tree 
and  a  tree  model,  which  embeds  the  prior  knowledge  about  the 
images.  Tree  matching  is  done  in  a  multilevel  way,  processing 
different  tree  levels  sequentially.  For  each  level,  an 
optimization  process  is  performed,  in  which  an  error  function, 
obtained  from  differences  between  the  model  and  the 
segmented  tree,  is  minimized.  13  parameters,  concerning  gray 
level,  shape,  position  and  connectivity,  are  used  to  characterize 
the  objects.  The  model  is  obtained  from  a  set  of  training 
images,  assigning  manual  labels  to  tree  nodes  with  a  user 
interface  designed  especially  for  this  purpose.  Three- 
dimensional,  multicomponent  images  can  be  processed  by 
adapting  gradient  and  parameter  calculation.  The  system  has 
been  tested  for  intracranial  cavity  segmentation  in  magnetic 
resonance  images,  giving  accurate  results. 

Keywords-  Medical  image  analysis,  model-driven  segmentation, 
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i.  Introduction 

Medical  image  segmentation  remains  an  open  problem 
in  many  different  applications.  Many  researchers  have  given 
specific,  ad  hoc  solutions  for  certain  problems.  General 
segmentation  schemes,  however,  are  far  from  being  accurate 
for  a  majority  of  cases  [1], 

It  is  unlikely  that  a  correct  automatic  segmentation 
could  be  obtained  without  relying  on  prior  knowledge  about 
the  images.  For  this  purpose,  the  utilization  of  anatomical 
models  is  a  good  way  to  represent  this  knowledge  [2],  [3], 

In  the  search  for  an  adequate  model  of  image  structure, 
multiscale  approaches  have  received  special  attention  [4], 
[5].  Compact  descriptions  can  be  achieved  by  using  singular 
points  of  the  image.  Several  research  studies  have  been 
carried  out  to  study  the  behaviour  of  these  points  under  the 
application  of  filters.  Local  extrema  of  the  image  can  be 
used  as  the  special  points  for  such  study:  combined  with 
mathematical  morphology  filters,  it  is  possible  to  obtain  a 
multiscale  description  [6],  providing  interesting  advantages: 
local  extrema  suffer  no  displacement  under  the  application 
of  dilation-erosion  filters;  furthermore,  a  complete 
segmentation  of  the  image,  at  different  scales,  is  easily 
obtained  by  watersheds  calculation  [7]. 

A  tree  description  is  especially  adequate  to  represent 
segmentations  at  different  scales.  This  trees  can  be  obtained 
efficiently  by  using  a  modification  of  the  immersion 
algorithm  for  watershed  calculation  [8].  Each  node  of  the 
tree  corresponds  to  a  minimum  of  the  gradient,  while  its 
height  in  the  tree  indicates  the  corresponding  scale. 


We  have  developed  a  system  that  obtains  a  complete 
segmentation  of  the  image  by  calculating  efficiently  its  tree 
representation,  matching  it  then  with  a  tree  model,  which 
contains  the  prior  knowledge  about  the  images.  All  the 
aspects  of  the  process,  including  model  calculation,  are 
explained  in  the  following  section. 

II.  Methodology 

A.  Previous  Steps 

Watersheds  are  calculated  on  the  gradient  of  the  original 
image.  For  multicomponent  three-dimensional  images,  the 
vector  gradient,  as  defined  in  [9]  is  used.  Calculation  of  the 
gradient  components  in  the  three  axes  is  carried  out  using 
the  Sobel  operator. 

The  original  images  are  filtered  using  an  anisotropic 
three-dimensional  filter  [10]  to  reduce  the  effect  of  noise. 
After  gradient  calculation,  a  dual  h-reconstruction  [6]  is 
applied  to  reduce  the  number  of  minima  of  the  gradient 
image,  which  corresponds  directly  to  the  number  of  regions 
that  result  from  the  watersheds  calculation. 

B.  Tree  Calculation 

The  immersion  algorithm  calculates  efficiently  the 
watersheds  of  the  image.  This  technique  uses  an  analogy 
with  the  gradual  immersion  of  a  topological  surface  (the 
image,  in  which  high  gray  levels  are  interpreted  as  high 
elevations,  and  vice  versa),  into  water.  We  have 
implemented  a  modification  of  this  technique,  in  which  the 
catchment  basins  merge  during  the  process  to  produce 
higher-level  regions,  as  shown  in  fig.  1  for  a  one¬ 
dimensional  signal  (two-  or  three-  dimensional  signals  are 
processed  in  the  same  way).  In  this  case,  the  normal 
immersion  process  would  have  produced  four  regions, 
labelled  as  1-4  in  fig.  1.  Our  modification  maintains  these 
regions,  while  adding  high-level  nodes  5  and  6.  To  avoid  an 
excessive  number  of  nodes,  new  ones  are  created  only  if 
necessary:  in  figure  1,  when  regions  5  and  6  merge,  a  link  is 
created  instead  of  a  new  node.  The  direction  of  this  link 
(that  is,  which  one  of  the  two  nodes  becomes  the  father) 
depends  on  the  extinction  value  of  the  corresponding 
minima,  which  in  fig.  1  corresponds  simply  to  the  size  of  the 
region  at  the  merging  time:  region  5  is  larger  in  this  case. 
Other  extinction  values,  such  as  the  dynamics  (depth  of  the 
minima)  can  also  be  applied  without  changing  the  tree 
building  algorithm. 
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Figure  1.  Construction  of  the  tree  during  the  immersion 
process 


Complete  segmentations  can  be  obtained  at  different 
scales,  by  taking  the  nodes  at  the  appropriate  level  and 
considering  the  regions  associated  with  each  one  of  them.  In 
fig.  1,  the  coarsest  scale  is  represented  only  by  node  5, 
which  covers  all  the  domain,  while  the  finest  one 
corresponds  to  the  output  of  the  unmodified  immersion 
procedure,  that  is,  nodes  1,  2,  3,  4. 

C.  Model  Obtaining 

Image  segmentation  is  achieved  by  transformation  of 
the  tree  obtained  by  the  modified  immersion  algorithm  into 
a  tree  model,  which  embeds  the  prior  knowledge  about 
anatomical  objects.  This  model  is  obtained  from  manual 
segmentation  of  a  set  of  training  images. 

Manual  segmentation  is  carried  out  in  the  following 
way:  a  tree,  corresponding  to  the  first  training  image  is 
constructed,  as  described.  Manual  labelling  of  the  obtained 
nodes  is  then  performed  using  an  interactive  user  interface. 
Labels  can  be  assigned  at  any  scale  level  of  the  tree,  thus 
minimizing  segmentation  time.  After  the  first  image  has 
been  accurately  segmented,  a  pruning  process  is  applied  to 
obtain  a  simplified  tree,  in  which  each  label  corresponds  to  a 
single  tree  node.  The  tree  topology  obtained  in  this  first 
segmentation  is  manually  validated  by  the  user,  and  can  be 
modified  to  eliminate  small  nodes,  which  may  appear  at 
undesired  locations  of  the  tree.  This  topology  is  preserved 
throughout  the  rest  of  the  training  process. 

Once  topology  has  been  determined,  the  following 
parameters  are  obtained  for  each  node  of  the  model  tree: 


•  Gray  level  parameters:  Mean  value,  standard  deviation 
and  maximum/minimum  values.  For  multicomponent 
images,  all  those  parameters  are  calculated  for  each 
image  component. 

•  Shape  parameters:  Area,  perimeter,  compactness, 
central  second-order  moments  p,20  and  po2,,  eccentricity 
and  orientation  of  the  major  axis  are  used.  In  case  we 
work  with  three-dimensional  images,  volume  and 
external  area  are  used  instead  of  area  and  perimeter, 
second-order  moments  used  are  p,2oo,  p.02o  and  Hoo2>  and 
orientations  of  the  main  and  the  second  axes  are  used. 

•  Node  centroid 

Parameters  are  calculated  for  all  the  training  images,  and 
their  mean  and  maximum  deviation  are  stored.  The  use  of 
these  values  is  explained  further  in  the  next  section. 

D.  Tree  Matching 

To  segment  a  new  image,  a  matching  process  between 
its  corresponding  tree,  calculated  as  described  in  section 
II. B,  and  the  model  for  this  type  of  images,  is  carried  out. 
The  usual  situation  is  the  following:  the  tree  obtained  from 
the  new  image,  which  contains  hundreds  of  nodes,  must  be 
transformed  into  a  much  simpler  model.  So  we  can  assume 
that  the  transformation  of  the  new  tree  into  the  model  is 
always  a  tree  simplification. 

This  simplification  problem  can  be  solved  by 
performing  a  classification:  the  nodes  of  the  tree  model  can 
be  seen  as  classes,  which  must  be  assigned  to  the  nodes  of 
the  new  tree.  Classification  is  carried  out  hierarchically, 
starting  at  the  coarsest  level  of  both  trees  and  descending 
step  by  step. 

At  each  level,  classification  is  carried  out  by  an 
optimisation  process,  in  which  an  error  function  is 
minimized.  This  error  function  quantifies  the  difference,  at 
the  current  scale,  between  the  tree  model  and  the 
classification,  and  is  defined  as  follows: 

Err  =  £  £ 3>(  |  fy  (new ) - ./.  (model)  | )  (1) 

f=l  2=1 

where  K  is  the  number  of  classes  (nodes  of  the  model  at  the 
current  scale),  F  is  the  total  number  of  parameters  used 
(listed  in  section  II. B),  fy  represents  the  value  of  parameter  j 
in  class  i,  and  <f>  is  the  weighting  function  that  is  shown  in 
fig.  2.  Two  parameters  must  be  determined:  Ay  |im,  which 
corresponds  to  the  maximum  “zero  error”  value,  and  the 
slope  of  the  function,  Wy.  They  are  determined  from  the 
manual  segmentations  of  the  training  images,  as  follows: 

•  fij(model)  is  the  mean  value  of  parameter  j  in  class  i  of 
the  images  in  the  training  set; 
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Figure  2.  Weighting  function  used  in 
error  calculation 


•  Ajjjim  receives  the  value  of  the  highest  difference  from 
the  mean  value  of  parameter  j  in  class  i .  In  this  way,  it 
is  assured  that  all  the  manual  segmentations  given  to 
the  training  images  have  zero  error. 

•  Wy  can  be  adjusted  to  assign  different  weights  to 
different  parameters  or  classes.  In  our  case,  wy  is 
assigned  the  inverse  of  Ay  |im,  so  that  differences 
produce  a  stronger  penalization  for  parameters  with  a 
smaller  spread. 

An  important  condition  that  should  be  achieved  is 
connectivity,  i.e.,  that  each  node  in  the  model  tree  should 
correspond  to  a  single  connected  component  of  the  labelled 
image.  So  we  define  the  connectivity  parameter  as  follows: 

A  ( i j 

^  /.s  max  comp  V1/ 

Con(i)  = - — —  (2) 

AO) 

where  Amax  comp(i)  is  the  area  of  the  maximum  connected 
component  of  node  i,  and  A(i)  the  total  area  of  this  node. 
For  three-dimensional  images,  area  is  substituted  by  volume 
in  (2). 

As  we  have  mentioned,  connectivity  should  always 
reach  a  value  of  1  at  the  end  of  the  optimisation  process.  A 
high  weight  of  the  connectivity  does  not  assure  that  this 
condition  is  always  reached  and,  moreover,  causes  local 
minima  of  the  error  function  to  appear.  To  prevent  these 
complications,  we  apply  initially  a  zero  weight  to  the 
connectivity,  until  the  process  converges  to  a  minimum. 
This  weight  is  then  progressively  incremented  until  total 
connectivity  is  achieved. 

The  optimization  process,  at  each  level,  consists  of  a 
series  of  movements  in  the  space  of  possible  classifications. 
Two  types  of  movements  are  allowed: 

•  Label  change:  assignation  of  a  different  label  to  a  node. 

•  Node  splitting:  one  of  the  nodes  is  partitioned  into 
different  classes,  by  taking  all  its  direct  descendants 
and  assigning  them  different  labels. 

Initial  classification  is  important,  as  it  can  affect  the 
process  in  two  ways:  avoiding  local  minima  of  the  error 


function  and  reducing  the  number  of  movements  necessary. 
Prior  knowledge  about  the  specific  segmentation  problem 
could  be  much  useful  to  determine  the  initialization  state. 
Anyway,  to  keep  the  segmentation  system  applicable  to  any 
type  of  medical  image,  we  take  as  initial  state  the  best  one  of 
a  number  of  random  solutions. 

The  optimization  technique  used  is  the  steepest  descent 
algorithm:  at  each  state,  all  its  neighbors  (that  is,  all  the 
steps  that  can  be  reached  with  a  single  movement)  are 
tested,  and  the  one  with  the  smallest  error  is  taken.  In  case 
none  of  the  neighbors  has  an  error  smaller  than  the  current 
value,  a  minimum  has  been  reached,  so  the  optimization 
process  stops. 

In  order  to  prevent  the  system  getting  stuck  in  a  local 
minimum,  several  initial  values,  taken  at  random  locations, 
are  tried.  The  final  solution  is  the  best  one  obtained  from 
these  different  initial  locations. 

m.  Results 

To  demonstrate  the  applicability  of  the  system  to 
medical  images,  it  has  been  used  to  segment  the  intracranial 
cavity  (ICC)  from  magnetic  resonance  images  of  patients 
with  multiple  sclerosis.  A  simple  tree  was  generated,  in 
which  the  head  is  first  extracted,  followed  by  a  classification 
in  three  classes:  skin,  bone  and  ICC.  The  tree  model  is 
shown  in  fig.  3. 

The  system  was  trained  using  10  images  of  different 
patients.  One  central  slice,  at  approximately  the  same 
location,  was  taken  from  each  study.  The  trained  tree  was 
then  used  to  segment  a  different  set  of  images.  Some  results 
are  shown  in  fig.  4,  where  head  and  ICC  contours  are 
superimposed  on  the  original  images:  as  can  be  seen,  results 
are  qualitatively  correct.  Cases  such  as  the  lower  left  of  fig. 
4,  where  cortical  bone  is  thin,  are  especially  challenging  to 
conventional,  gray-based  only  techniques,  as  it  is  difficult  to 
separate  ICC  from  bone.  In  our  system,  the  introduction  of 
shape  parameters  and  the  minimization  of  a  global  function 
allow  to  obtain  correct  results  even  in  these  cases. 
Calculation  time  needed  to  process  the  images  in  fig.  4  was 
between  1  and  3  minutes  in  a  PC  Pentium  III. 


Figure  3.  Tree  model  used  for 
segmentation  of  intracranial  cavity 
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Figure  4.  Results  for  head  and  ICC  contours  of  different 
patients 

iv.  Discussion 


Many  different  applications  of  medical  images  require 
segmentation  of  the  anatomical  structures.  Specific  systems 
have  been  developed  for  special  applications.  However,  if 
segmentation  of  a  different  object  is  needed,  or  even  if  the 
acquisition  parameters  vary,  application-specific  systems 
must  be  deeply  redesigned  for  the  new  application 
We  have  developed  a  system  that  separates  the  segmentation 
process  from  the  a  priori  knowledge  of  the  images  (the 
model),  so  that,  if  we  have  to  segment  a  different  type  of 
image,  only  this  model  must  be  recalculated.  Furthermore, 
we  have  designed  an  efficient  way  to  determine  the  models, 
starting  from  the  original  tree.  This  model  calculation 
interface  can  also  be  used  independently  as  a  semiautomatic 
segmentation  system. 

The  system  can  work  on  three-dimensional, 
multicomponent  images.  Processing  time,  however,  may  be 
a  problem  if  very  big  data  volumes  are  used.  Alternatives  for 
parallel  calculation  of  the  data  are  being  studied. 

The  system  presented  has  demonstrated  its  performance 
in  segmentation  of  the  intracranial  cavity  from  MR  images. 
A  detailed  quantitative  evaluation  is  being  carried  out  at  the 
moment.  Different  applications,  such  as  liver  segmentation 
from  CT-MR  images,  will  also  be  tested. 


iv.  Conclusion 

A  new,  general  purpose  medical  image  segmentation 
system  has  been  presented.  The  use  of  a  tree  model  makes 
the  system  easily  applicable  to  different  anatomical 
locations  or  modalities.  The  system  works  on  3-D, 
multicomponent  images,  thus  taking  profit  of  multimodal 
patient  explorations.  Initial  results  have  demonstrated  the 
applicability  of  the  system  to  specific  medical  problems. 
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